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SUMMARY 


The motion of an elastically supported cylinder forced by a nonlinear, quasi-static, aerodynamic 
model with the unusual feature of a motion-dependent forcing frequency has been studied. Numerical 
solutions for the motion and the Lyapunov exponents are presented for three forcing amplitudes and two 
frequencies (1.0 and 1.1 times the Strouhal frequency). Initially, positive Lyapunov exponents occur and 
the motion can appear chaotic. After thousands of characteristic times, the motion changes to a motion 
(verified analytically) that is periodic and damped. This periodic, damped motion has not been observed 
experimentally, thus raising questions concerning the modeling. 

INTRODUCTION 


Vortex excitation of elastic structures such as aircraft wings, helicopter blades, and cables can lead 
to a wide diversity of motions. Experiments involving simple configurations, that of cylinders or wires 
in air or water, demonstrate the range and complexity of motions yielded by the nonlinear coupling of 
mechanical and fluid motions. Sreenivasan (ref. 1) and Van Atta and Gharib (ref. 2) showed that, for 
cylinders suspended in fluid, regions of periodic, quasi-periodic, and chaotic motion exist. In addition 
their experiments, as well as those reported in a review article by Bearman (ref. 3), indicate that for certain 
Reynolds numbers the shedding frequency locks onto the natural frequency of the cylinder or onto one of 
its subharmonics. Feng’s experiments (ref. 4) show that, when the shedding frequency and the frequency of 
the motion of the cylinder are locked together, the amplitude of the motion is large and exhibits hysteresis. 

In this paper we investigate the ability of the one-degree-of-ffeedom, nonlinear, quasi-static model 
proposed by Tobak et al. (ref. 5) to describe the range of motions observed experimentally for an elasti- 
cally supported cylinder forced by vortex shedding. In this model, the mechanical and fluid motions are 
coupled through a motion-dependent forcing frequency that models the dependency of the vortex shedding 
frequency on the motion of the cylinder. This motion-dependent forcing frequency is unique in terms of 
previous studies of nonlinear systems. 

This one-degree-of-freedom model with its unusual forcing frequency produces some unexpected 
and unusual results. Solutions that for thousands or even millions of characteristic times may appear peri- 
odic, quasi-periodic, or chaotic eventually settle into a long-term motion that is damped and periodic. This 
long-term motion has not been observed experimentally and does not appear to be physically reasonable. 
Once the motion has assumed its long-term behavior, the computed Lyapunov exponents are unbounded. 


1 


Both the damped, periodic, long-term motion and the unbounded nature of the Lyapunov exponents are 
verified analytically and result directly from the motion-dependent forcing frequency. 

EQUATION OF MOTION 


Since the motion of an elastically supported cylinder in a fluid with a constant free-stream velocity 
Uo is primarily in a direction perpendicular to the free stream (Bearman, ref. 3), Tobak et al. (ref. 5) used 
a one-degree-of-freedom system for modeling the motion: 

Mh, + 2ph + k 2 h = N\ + N 2 

Here, M is the mass per unit length of the cylinder, p is the spring damping constant, n 2 is the spring 
constant, h is the position of the cylinder perpendicular to the free stream, and N\ and N 2 are aerodynamic 
forces acting upon the cylinder in a direction normal to the free stream (fig. 1). Furthermore, they assumed 
that, in the Reynolds number ranges where vortices are shed alternately from each side of the cylinder, the 
pressure fluctuations at the surface of the cylinder are proportional to sin (u/t) where oj is the shedding 
frequency and that the pressure can be written as the Fourier series 

^\„(r) cos (nff) + (r) sin (nB) sin (wt) 

_ n n 


Here, Uoo is the free-stream velocity observed at the cylinder and is equal to yU^ + h 2 . 


Pa ~ 2 PooUqqA 


Using this expression for the pressure and keeping only first-order terms gives the following expres- 
sions for the drag per unit length of the cylinder, D , and lift per unit length, L : 

D = jPooUlDCo 
and 

L = ^PooUIDCl sin (ut) 

The drag and lift coefficients, Cd and Cl, are assumed constant and do not depend upon the Reynolds 
number. The components of drag and lift in the direction perpendicular to the free stream are 


and 


N\ = D sin a 


N 2 = L cos a 


where ex = tan -1 (h/Uo) is the local angle of attack. Using these expressions for N \ , N 2 , and a gives the 
following dimensional form of the equation of motion: 

Mh +2 ph+ K 2 h = —pooUlcD 


Cdjj— + Cl - jrj— sin (ut) 

U 00 U 00 
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The nondimensional equation of motion is obtained by letting h — > Dh andt — > tD/Uo, giving 
2 L %Dh Uo . ( D \ 

Cd ~u^* CL uZ sm Vn t ) 


Dh + 2^Dh + D^h = lp m UlD 


or 


l . „ D L . 1 ( kD\ 2 1 Poo D 2 U Uoo L . „ Uoo . ( D ^ W 

h + 2, W + m { «T ) h = T — [° D tfc' k ' * Cl W “M 


where 


f/oo = \/c/ 0 2 + [^o/^p/i] 2 
= £7 0 >/l + T 2 

Nondimensionalizing p and « by — > p / and — * k', respectively, gives 


h + 2 p'A + k' 2 h = eV^ 1 + A 2 


/ D N 
Cz?/i + Cl sin w — t 
\ Uo , 


where e = \p (x D 1 /M is proportional to the ratio of the density of the fluid and the density of the cylinder. 
The dimensional vortex shedding frequency, w, is assumed by Tobak et al. (ref. 5) to be proportional to 
k a , the Strouhal number 

w = k s Uoo/D 

The Strouhal number is independent of the Reynolds number and assumed to be constant, k a = 0 An. 
Using this expression for to and dropping the prime notation gives the nondimensional form of the equation 
of motion: 

h + 2p/i + K? h = c\/ 1 + ^C/j/i + Cl sin ^ k 3 tV 1 + /i 2 )] 


The forcing frequency in this last equation is velocity dependent, and the forcing term cannot be 
written as a simple Fourier time series. An unusual and perhaps unrealistic feature of this motion-dependent 
forcing frequency is that unless h is constant, a shift in the time origin will change the forcing function. 

To demonstrate this point, sin (k a t vTThT^ has been plotted in figure 2 for two different time intervals, 

0 < t < An/k a and 1000 < t < 1000+ 4ir/k 3 ] h has been assumed to be sinusoidal, h = 0.1 sin (k a t). 
For small times, the dependency of the forcing term upon the motion can be ignored, and the forcing term 
appears sinusoidal. As time increases, this is no longer the case. 

The effect of time on the forcing function can be explained by approximating the forcing argument 
by k a t + 0 .005 k a t sin 2 ( k a t) . For small times, the second term in the argument can be ignored; it is small 
with respect to the first term and with respect to 2 n. As time increases, the second term remains small 
with respect to the first term, but it varies from 0 to a value that is larger than, and can be much larger than, 
2 7r. This causes the forcing function to pass through 0 several times in the interval to < t < to + 2n/k a . 
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The number of zeroes in the interval to < t < to + 2 rr/fc, will increase as to increases, and, therefore, the 
forcing function has a dependency upon the time origin. This dependency of the forcing function upon the 
time origin will also effect the time history of the motion, as will be described in later sections. 

NUMERICAL INTEGRATION 


The nonlinear equation derived above can be written as a third-order autonomous system. Let 

/ h \ 

h 


y = 


\k s tVTT¥/ 


Then 


dy 

dt 


/ 


h 

h 


\ 


\ k a V 1 + /i 2 + k a thh J 
V VI + h 2/ 

or, in terms of the y, ’s, the third-order autonomous system is 


dy_ 

dt 


( 


J/2 


\ 


-2^1/2 - « 2 yi + €y/l + y2 2 [CnV 2 + Cl sin(y 3 )] 

\k 3 \/\ + y 2 2 + |-2 /iy 2 - n 2 y\ + e\/l + y2 2 [CdVi + C L sin(y 3 )]| ) 

y 2 

The Jacobian of this system is 


/ 0 


0 


\ 


-K 2 022 e Ci\fi + 1/2 2 cos(y 3 ) 




0 3 2 


033 


/ 


where 022^032, and a 33 are given by the following complicated expressions: 


022 = — 2 /J + 


/ [Cd (l + 2y 2 ) + Cz,y2 sin ( y 3 ) ] 

V 1 + 


Q 3 2 — 


1/2 


\/l + yr 


-2/i + 


vl + yr 


[C D (1 + 2y 2 ) + C L y 2 sin(y 3 )] J 


+ {-2/XJ/2 - ^ 2 yi + €\/l + y2 2 [ Czm + Cx sin(y 3 )] j YJ 


1/3(1 -yf) 


+ y 2 z ) 


2 A 2 


and 
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1/2 


033 = |-2^y 2 - « 2 yi + fV 1 + yi 2 [CdV 2 + Cl sin(y 3 )] j y^- 


Vi 


e\/r + 1/2 2 Cl cos (y 3 ) 


j/2 i/3 

i + y\ 


This third-order autonomous system of equations was integrated using the Bulirsch-Stoer numerical 
estimator (ref. 6). This integration scheme divides the integration interval, At = U+\ — f,, into an even num- 
ber of intervals and calculates a value for y( ti+\ ) by using modified midpoint estimates at t ,+ M A t/( 2 N ) , 
M = 1,2 N. Each modified midpoint estimate is second-order accurate in [ A t/( 2 N) ] , and the calculated 
value of y(ti+\ ) , therefore, is dependent upon the number of intervals: 

y( t i+l ) =>Y(t i+u At/(2N)) 

Rational extrapolation of the function F(f,+i ,At/(2N)) to Y ( 1 , 0) gives an estimate for y(t, +1 ) us- 
ing an infinite number of midpoints and, hence, an infinitesimal stepsize. An error is estimated for the 
extrapolated value of y{ U +\ ) , and if this estimated error does not fall within given bounds, N is increased. 
This integration scheme is typically more accurate than higher-order Runge-Kutta methods and allows a 
built in limit on the estimated error size. 


The Lyapunov exponents were calculated using the method of Wolf et al. (ref. 7) and the Bulirsch- 
Stoer integration scheme. In Wolf’s method, an orthonormal set of vectors, x ; (f,), j = 1 , 3 is defined at 
t{. At fi+i , these vectors are transformed as follows: 
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Xj(ti+l) = Xj(t{) + 



Jxj(t)dt 


where J is the Jacobian of the autonomous set of equations. The values of the Lyapunov exponents at t, 
are found by equating the length of xi(£,+i) to exp[ Ai.dfj+i — t,)], the area spanned by xi(t,+i) and 
X 2 ( t»+i ) to exp [(>]_,• + A 2) ,-)( t i+i - U ) ] , and the volume spanned by x\ ( t,+i ) , x 2 (f,+i ) , and x 3 ( f,+i ) to 
exp [ ( \ i t i + A 2 + A 3 ,» ) ( t i + 1 — t i) ] . The vectors x ; - ( t ,+ 1 ) are then orthonormalized using the Gram-Schmidt 
procedure, and the whole process is repeated. The Lyapunov exponents \j are defined by 




n — 1 /*^» 

^ — ti) / 

«'= 0 J ' to 


( t)dt 


tn to 


tn to 


Note that the Lyapunov exponents as defined here will approach a constant value only if the A ; /s approach 
a constant value or oscillate about a constant mean value as t, — > oo. If, on the other hand, X ;> ,- is time 
dependent and of the form = A + Bt f , the previous expression gives 


\j = A + 


b 

1 + 13 t n -to 


5 



or, for t n > to , 


\j = A + 


B 

1 + 0 


tS 


Although the exponent of t and the sign of the term are properly captured in the above expression, the 
coefficient is not. A multiplicative correction to the computed Lyapunov exponents must be considered if 
exact values are required. 


RESULTS 


Three sets of values were chosen for the nondimensional damping coefficient /i, the aerodynamic 
damping parameter eCp, and the aerodynamic forcing parameter eCx • For the large-forcing cases, //, cCn, 
and eCx were chosen as 0 .0054 , — 1 , and +1 , respectively. The natural frequencies were k = 0 .2033 k s , 
k = k s , and k = 1.1 k a . Although these cases did not correspond to any available experimental data, they 
were used because their transients are short-lived and the long-term motion is rapidly assumed. A large- 
forcing case with small damping and a large natural frequency (eCx = 1 , eCo = -0 .001 , /i = 0 .001 , and 
k = 2 .3033 k s ) was used to further investigate the mechanisms of transition from the short-term, transient 
behavior to the long-term solution. 

For the medium-forcing cases, the parameters were chosen to approximate Feng’s (ref. 4) data for a 
circular cylinder in air and are given by /i = 0 .0054 , cCd = —0 .0294 , and eCx = 0 .0147 . The natural 
frequencies were k = 0 .2033 k 3 , k = k s , and k = 1 . 1 k a . For the small-forcing cases, /i, eCn, and cGl 
were chosen as 0.00001, -0.00020, and 0.00007, respectively, to approximate the values for Van Atta and 
Gharib’s experiments using piano wires in air (ref. 2). The natural frequencies were k = k 3 and n = 1.1 k a . 

The computations were all performed using a CRAY XMR Millions of integration steps were nec- 
essary for the study of the short-term motion, the long-term motion, and the transition from the short-term 
to the long-term motion. The integration step size depended greatly upon the forcing amplitude and only 
slightly upon the natural frequency of the spring and time. The value of the step size was automatically var- 
ied between a user-imposed upper and lower bound as the integration proceeded so that rapid oscillations 
and sudden jumps in the acceleration and forcing function could be accurately computed. The lower bound 
on the step size was arbitrarily chosen to be 10 -6 (the CRAY XMP has 15-digit accuracy; hence, this step 
size would still be reasonable for times of order 10 6 ). The upper bound varied from a maximum of 7t/80 
(approximately 125 points/cycle) for the small-forcing cases to tt/ 640 (approximately 1000 points/cycle) 
for the large-forcing cases. 


Large Forcing 

For three of the four large-forcing cases, those with natural frequencies k = 0 .2033 k s , k = k s , and 
k = 1 . 1 k a , the motion rapidly assumed a long-term behavior that was independent of the natural frequency 
and was periodic and damped. One of these three cases, the « = 0 .2033 k s case, is used here to demonstrate 
the features of the long-term behavior. 
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In figure 3, plots of the long-term motion versus time, the phase portraits (h versus h), and the peak 
amplitude versus time are shown for the k = 0 .2033 case. The plots of the motion versus time and the 
phase portraits contain similar information; the amplitude of the motion, the frequency of the motion, and 
the rapidity of changes in the motion can be determined from either the plots of motion versus time or 
the phase portraits. However, the phase of the motion cannot be determined from the phase portraits, and 
determining the frequency of the motion and the rapidity of changes in the motion involves an integration 
along the ( h, h) curve. Therefore, in this paper, the discussion of the motion will center on the motion 
versus time plots and not on the phase portraits. 

In figure 3a, plots of the long-term position h, velocity h, and the acceleration h indicate that the 
frequency of the motion is equal to the Strouhal frequency, as would be predicted by a linearized form 
of the equation of motion obtained by approximating the forcing frequency by its average value and by 
assuming that h < Cl- However, the amplitudes of h and h are much smaller than that predicted by 
the linear equation, and, furthermore, the motion is damped and appears to be approaching a fixed point, 

(M) = (0,0). 

This appears to contradict the fact that, for finite times, no fixed point solution exists. Letting h, h, 
and h go to zero in the equation of motion gives 

0 = cCl sin(k 3 t) 

This equation can not be satisfied for all t unless k 3 = 0 , and, therefore, no fixed point solution exists for 
k 3 y 0 . However, consider another approach to this problem. Assume that, for large times, the motion is 
such that the forcing function is minimized. This implies that 

sin (k 3 t\/ 1 + A 2 j = 0 

or 

k 3 t\/l + h 2 = rm 

giving 



where n is an integer. This expression is valid only for t < rm/k s . For t > rm/ k s , n must be replaced by 
a larger integer, n —> n+ 1 . Iterating this procedure gives 




(n+ 1) 7r 

k 3 t 


2 


- 1 


for 


rm , ^ 

feT <t - 


(» + l)7r 

k 3 


By assuming a sign change when n — > n+ 1 (this is supported by the computations), the general form for 
the velocity at large times can be written as 


h 


±( 



(m — 1) 7r raw 

for — k <t -lT’ 

ft's ft's 


m = no , oo 
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where no is proportional to the time at which the long-term behavior began. In figure 4, the previous 
function is plotted for one of the time ranges shown in figure 3a. Note that the calculated velocity in 
figure 3a and the conjectured velocity in figure 4 show the same characteristics; the period is 2n/k s and 
there are jumps whenever h approaches 0. The amplitudes of the calculated and conjectured velocity are 
virtually identical. Any discrepancies between the calculated velocity and conjectured velocity may be due 
to a variety of reasons: 1) the conjectured velocity only solves the equation of motion exactly as t — ► oo; 
2) for the calculated motion there may be a small offset in the time origin (t — > t — to ) that decreases in 
importance as time increases; 3) the argument of the calculated forcing function may have a small departure 
from rmr, k a tv 1 + fi 2 = run + 6/(f); or 4) the jumps occur over a finite time interval for the calculated 
velocity and an infinitesimal time interval for the conjectured velocity. 


For the conjectured velocity, the height of the jumps can be analytically determined as follows. The 
jump occurs when |/i| — > 0 and the argument of the forcing function changes from ( m — 1) 7r to m'n. Both 
the jump in the velocity and the amplitude of the velocity are given by lim^o h((m — 1) 7r//c s + |S|), or 



where m has been assumed to be large. The integer m is proportional to time and can be approximated by 
k 3 t/ 7 r. Therefore, both the velocity jump and amplitude can be written in terms of time as 


h 


max ^ 



2.24r 1/2 


and h converges to 0 at a rate f -1 / 2 . 


When jumps in the velocity occur, the cylinder is either at minimum or maximum displacement. 
Hence, it is possible to calculate the total displacement, A h = h max — h min , during the interval (m — 
\)Tr/k s <t < mix /k a by integrating the conjectured velocity; 


This gives 


A h = 


L 


Tmr/k s 

m-l)n/k s 


hdt 


l'm'n/k i 
J (m-\)n/k 3 



- 1 dt 



mix , 

+ 5X ln 
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Letting e = -^1 — ( rr ^-) 2 and expanding in terms of e gives 


A h = 


or 


7B7T£ 

ran 

k s 


2k s 

7717T 
+ 

2k a 

( E - 

7Tl7r i 

it' 

( £ 3 
^3 

+ . . . 


f 


7717T 1 

1 

1 - 

?r 

Ce 

A 

3 





e 

~ B ~Y 


e 2 e 3 e 4 


3 4 


m — 1 


TO 


-.3/2 


Assuming that m is large and, therefore, £< 1, the displacement is given by 

A , 2tt fl 

Ahp * TTV~ 

3k, V ro 


2 

Jk 


7T FTt\ _ 

V 


3 12r x ) 2 


Since \h m jn | = |/i ma x|, the amplitude of the conjectured position is given by 1 .86 £ -1 / 2 . As is the case 
with the conjectured velocity, the amplitude of h is proportional to f -1 / 2 , and h converges to 0. These 
conjectured results for the amplitudes are in good agreement with those computed for the large-forcing, 
low natural-frequency case (fig. 3c). 


The acceleration h does not go to a fixed point. As h and h go to zero, any changes in the forcing 
function must be canceled out by the acceleration. Hence, when the argument of the forcing function jumps 
from tott to ( to + 1) 7r and the forcing function jumps from 0 to ±cCl and back to 0 , h must also jump 
from 0 to ±eCz, and back to 0 . The derivative of the conjectured velocity for large t is 

h = k s (^) 1 

\ M / \/(TO7r) 2 — (A ; a t) 2 

pa — = 

Try to 2 — (to — 1 + a) 2 


where k s t = ( to — 1 + a) 7r, 0 < a < 1 . For large to. 


7r\/2m\/l — a 

For any given a 1 , the conjectured acceleration goes to zero, as do all other terms in the nonlinear 
equation of motion. At a = 1 , jumps occur in the velocity and the acceleration is undefined. However, if 
it is assumed that the jump in the velocity occurs over a finite time increment, that the magnitude of the 
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jump is given by k s t ) , and that the maximum acceleration is approximately equal to the amplitude 

of the forcing function, eCi, then the time interval over which the jump occurs, Af, can be approximated 
by 


Af K, 


1 /2 7 T 

eCz, V k a t 


This interval also decreases as . In figure 3 a, the plots of the computed acceleration do show that the 
acceleration changes from 0 to ±eCz, and back to 0 in a very short time interval and that, as time increases, 
the length of this interval decreases. 


Although these conjectured forms of h, h, and h, do not solve the equation of motion exactly for 
finite times, they do give insight as to the behavior and damping of the long-term solution. Instead of a 
motion with a nondecaying amplitude that forces higher frequencies into the motion as time increases, the 
long-term motion succeeds in minimizing the forcing function. Because of this minimization, the forcing 
function and acceleration are zero, except during very short time intervals which decrease in length as time 
increases. Since the peak forcing amplitude and peak acceleration are bounded and the length of the time 
interval during which the forcing function and acceleration are nonzero goes to zero, the acceleration and 
the forcing function do not impart any energy to the cylinder, and the motion dies out. 


For any choice of parameters (as will be discussed in the following sections), the motion eventually 
assumes the damped, periodic behavior just described. However, there has been no experimental evidence 
corroborating this damped, long-term motion. The experimental motions do not gradually die off, they 
do not show periodic jumps in the velocity, and they do not show a minimization of the acceleration. 
Therefore, this conjectured (and computed) long-term motion is nonphysical and will be referred to as 
such. 


This interpretation of the long-term behavior ( h and h fixed, j/3 = k 3 t varying) would seem to indicate 
that the motion should alternate between two fixed points; h = 0 , h = 0 , and j/3 mod2 n = 0 and h = 0 , 
h = 0 , and 1/3 mod2 7t = 7t. The solution periodically jumps from one fixed point to the other, and the 
period of the motion is given by the Strouhal frequency. This implies that one of the Lyapunov exponents 
should tend to zero; the other two should be negative. 

In figure 5 the Lyapunov exponents and the sum of the Lyapunov exponents are plotted for n = 
0 .2033 k a . Note that two of the exponents are small and negative and in the limit appear to approach zero. 
This contradicts the previous conjecture; however, these two Lyapunov exponents do remain negative for 
finite times, implying that the phase space does continue to contract. They do not reach zero, and h and h do 
not become fixed, for finite times. The fact that two of the Lyapunov exponents appear to be approaching 
zero as opposed to only one of the exponents may be a consequence of the slow convergence rate as well 
as the discontinuity in the asymptotic solution; it is not because the long-term solution actually falls on a 
two-torus. 

The third Lyapunov exponent and the sum of the Lyapunov exponents is negative and unbounded. 
The sum of the Lyapunov exponents represents the growth rate of the volume of the phase space, 
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V(t) oc exp (53 A,t). The volume at to + dt can also be written in terms of the volume at to and the 
trace of the Jacobian: 


V(t 0 + dt) = (\ + Tr Jdt)V(to) 
= exp (TV Jdf) V(fo) 


Iterating this equation and taking the limit dt — > 0 gives 

rti 


y(t 1 )=exp (J TrJdtjV(to) 
and the sum of the Lyapunov exponents is given by the expression 

V(U) 


V(t o) 


= exp 


- to) 

L«=i 


(£ 


= exp / TV J dt 


or 


A /‘'TV Jett 

i=l 


tl — to 

For large t (large 1 / 3 ), the trace of the Jacobian can be approximated by 

TV J ^ e\/l + y 2 2 C L cos (y 3 ) 

1 + y 2 

Using the long-term approximation for h, letting to = rmr/k a and t\ = rm/k a , and noting that cos(y 3 ) = 
cos( j'n) = ( — 1 ) ; gives 


]=m 


J 7T 


- 1 dt 


For n >> m 1 


= ±^ £ CW2^ 

1=1 

= ±^-cCL\/27rA:ati = ±1.25t ,//2 

and the sum of the Lyapunov exponents should be proportional to t 1 / 2 . Figure 5 shows that the corrected 
sum of the Lyapunov exponents is equal to 1 ,4t° 52 ; the deviation from t° 5 , as well as the error in the 
coefficient tCi^y/luka, may result from the use of the long-term (t — > 00 ) approximation over finite 
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times. As the time increases, the terms that are neglected for small times should become less important and 
the conjectured long-term behavior will better describe the motion and the sum of the Lyapunov exponents. 

Large Natural Frequency 

The large-forcing case with a large natural frequency and weak damping also went to the long-term 
behavior described previously. However, much more time was required for transition to the long-term 
behavior, and the transients showed an unusual behavior. In figure 6a, a time series for this case is shown. 
The amplitudes of the velocity and position are greater than those predicted by a linear system with damping 

and forcing coefficients given by 2/r + \zCd\ ( 1 + h 2 J and eCx ( 1 + h, 2 J , respectively, and they 
are much larger than the amplitudes predicted for the long-term behavior. 

The large velocity amplitude causes the argument of the forcing function through many multiples 
of 2 tt. For instance, consider the motion at f = 1550. The ampli tude of the velocity is 0.7861, and the 
argument of the forcing term can take on values ranging from k 3 tV\ + h 2 & 1950 for h - 0 to 2480. In 
the time that it takes the velocity to go from 0 to its peak amplitude (A t « 0 .5 ), the argument of the forcing 
function will change by approximately 1687T and the forcing function will go through zero approximately 
168 times. Since the oscillations caused by the variation in h, are rapid, they have no appreciable effect 
upon the velocity or position, and the velocity and position appear to be well-behaved periodic functions 
of time. 

The amplitudes of the velocity and position distributions decrease with time, but at different rates 
(fig. 6b). This difference in the damping rate is due to a change in the frequency of the motion and the 
departure of the motion from a simple sinusoidal motion. The initial motion of the cylinder has a frequency 
close to that of the natural frequency; as time increases, the frequency slowly decreases until it drops below 
the Strouhal frequency. This frequency variation is not given by the average value or maximum value of 
k 3 V\ + A . 2 ; a more complicated mechanism appears to be responsible for the frequency variations. 

As the frequency of the motion continues to decrease, the amplitude of the motion drastically de- 
creases and the long-term motion takes over. All the characteristics of the long-term motion, such as 
minimization of the acceleration, jumps in the velocity, and a very slow damping rate, are observed. There 
is what appears to be a small sinusoidal term superimposed upon the minimized forcing function. This 
may be due to some residuals from the short-term motion, and may disappear as the time increases. 

Even before the motion has settled into the long-term behavior, there appears to be some minimiza- 
tion of the acceleration. Obviously, in regions where the forcing function is oscillating rapidly, the ac- 
celeration must also oscillate rapidly; no other terms in the equation of motion can cancel out these rapid 
oscillations. However, during each cycle there are regions where the velocity is approximately constant 
and the acceleration is equal to zero. The length of these regions increases as time increases and as the fre- 
quency of the motion decreases. Furthermore, as time increases, the length of the high oscillation regions 
decreases, and the number of oscillations also decreases. Eventually, the high-frequency oscillations die 
out, and the acceleration and forcing function take on the characteristics of the long-term motion. 


12 



When the motion finally does assume the long-term behavior, there is a sudden change in the behavior 
of the Lyapunov exponents (fig. 6c); the small, monotonically increasing, positive Lyapunov exponent 
begins to decrease, and the small negative Lyapunov exponent becomes more negative. It appears that 
the sum of these two exponents will eventually become negative, and, possibly, both of these coefficients 
will become negative. Therefore, the long-term solution, as predicted by the Lyapunov exponents, lies 
at most on a one-torus, and possibly is a fixed point. The large negative Lyapunov exponent becomes 
less negative, but remains large. It is expected that once the contribution of the transient solution to the 
Lyapunov exponents becomes less pronounced, the sum of the Lyapunov exponents will assume the t 1 / 2 
behavior. 

Numerical Error and Long Time Series 

The previous computations for the large-forcing, large natural-frequency case indicate that at least 
two types of solutions exist for this equation of motion. The first is a sinusoidal, long-lived transient 
with a time-varying amplitude and frequency. The second is the nonphysical long-term motion with a 
minimized forcing function and acceleration. The number of time steps necessary for transition from the 
sinusoidal motion to the long-term motion is in the tens of millions, and the cumulative error may not be 
small. The amplitude, phase, and even character of the motion may be affected by this cumulative error, 
and the transition to the long-term motion may result from this error. The calculated motion may bear 
no resemblance to the exact solution for any given set of initial conditions. Furthermore, the cumulative 
error may also affect the sizes, signs, and trends of the Lyapunov exponents making them useless in the 
interpretation of the motion. 

Two problems in particular arise when using numerical solutions for chaotic systems. The first in- 
volves the dependency of the solution behavior on the time step. Tongue (ref. 8) showed for the Duffing 
equation that too large a step size could produce chaotic solutions instead of the correct periodic solu- 
tions. The second problem connected with numerical solutions of chaotic systems is the divergence of 
two solutions started arbitrarily close together. These solutions exponentially diverge in time (separation 
oc e Xit , X, ’s are the Lyapunov exponents), and two solutions that are similar for short periods of time may 
be vastly different after extended periods of time. For any numerical solution, small errors are introduced 
at each time step. Hence, for a chaotic system, the numerical solution may diverge drastically from the 
exact solution, and the numerical solution may even approach the wrong attractor. Therefore, for any nu- 
merical study of a chaotic system, the effects of both the step size and small errors on the solution must be 
investigated. 

Two methods were used to verify the computational solutions. The first involved changing the step 
size and the number of significant digits used in the computations. The results for the large-forcing, large 
natural-frequency case using two different user-defined maximum step sizes, dt = tt/1600 and dt = 
tt/6400, are shown in figure 7. The double-precision results for dt = tt/1600 are also shown in figure 7. 
All three calculations were started at t = 6991 with h = 0 .25 and h = 0 . 

The amplitude, period, and phase of the motion do not appear to be dependent upon the maximum 
allowed step size or the number of significant digits used in the computation. Also, the qualitative features 
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of the computed Lyapunov exponents is unaffected by the change in the maximum allowed step size or by 
the change from single to double precision. For all three cases, the first Lyapunov exponent is small, posi- 
tive, and has regions of large frequency oscillations, the second Lyapunov exponent is small and negative 
with periodic spikes, and the third Lyapunov exponent is large, negative, and composed of a periodic term 
superimposed upon a monotonically decreasing term. 

Although the Lyapunov exponents for these three cases do share the same general features, there are 
some quantitative differences. For the third Lyapunov exponent, the slope of the monotonically decreasing 
term is dependent upon the step size, indicating that there is some error in the damping rate of the motion. 
The minimum values of the second exponent as well as the height of the spikes also appear to be dependent 
upon the step size. The quantitative differences between the single- and double-precision solutions were 
minimal. 

It should be noted that these values for the Lyapunov exponents are vastly different from those of 
the solution started at t = 0 . Since the calculation of the Lyapunov exponents involves a summation 
over time, this is not unexpected; the values of the Lyapunov exponents for two solutions started at two 
different times would only agree if the exponents are constant or periodic for almost the entire time span 
of the calculations. 

The second method of verifying the calculated motion involved checking the effect of amplitude 
deviations on the general behavior of the solution. Solutions for the large-forcing, small-damping case 
were started at three times, t = 99 l,f = 6,991, and t = 13,991 with three different starting amplitudes, 
h = 0 .0 , h = 0 . 1 , and h = 0 .2 ; the initial value for h was chosen to be zero. These amplitudes were of 
the order of, but not identical to, the amplitudes calculated for the solution started at t = 0 , and the results 
of these computations are shown in figures 8a and 8b. 

When the starting amplitude was large, h = 0.2, the motion took on the behavior of the short term, 
periodic motion. For no initial displacement or velocity, the long-term motion was immediately assumed, 
and the amplitude of the long-term motion was approximately equal to that predicted by analysis. For 
the initial displacement equal to 0.1, the motion exhibited all the characteristics of the long-term motion; 
the forcing function was minimized, the frequency of the motion was equal to the Strouhal frequency, 
the velocity plots showed sudden jumps, and the amplitude of the velocity was approximately equal to 
that predicted by analysis. Added to the conjectured long-term position was a monotonically decreasing 
function which allowed the initial conditions to be satisfied. 

In figure 9, the motion is plotted for an initial time equal to 6991 and three initial displacements, 
h = 0.15, 0.20, and 0.25. Note that the computed motions for these three cases show all the short- 
term characteristics: the motion is periodic, there are regions of constant velocity and minimization of 
the acceleration, and there are regions where large frequency oscillations occur in the acceleration. The 
amplitude and phase of the motion is dependent upon the initial conditions, and the frequency of the motion 
is primarily amplitude dependent. 
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The results shown in figures 8 and 9 indicate that there are at least two types of solutions, the periodic 
long-term solution with a frequency equal to the Strouhal frequency and a periodic short-term solution 
with a frequency that is dependent upon the amplitude of the motion. The general characteristics of these 
two solutions do not change when there are small errors in amplitude or displacement. The character of 
the solution does change when the motion transforms from the short-term to long-term behavior. This 
transition depends largely upon the amplitude of the motion and only slightly upon the time; therefore, 
errors in the amplitude could induce the transition. However, since the forcing function can not sustain 
a periodic motion, it can be argued that the motion must slowly damp out. Hence, at some point, the 
amplitude of the motion should be small enough that the motion will assume the long-term behavior. The 
time that this transition occurs may be dependent upon step size and the accuracy of the solution, but the 
transition will occur. 

These two methods of solution verification (i.e., the study of the dependence of the solution on the 
step size and number of significant digits and the investigation of the effect of small amplitude changes on 
the solution) seem to indicate that the computed solutions mimic the behavior of the exact solution. Al- 
though these solutions may not give exact amplitudes, frequencies, times of transition to the long-term 
behavior, or Lyapunov exponents, they can be used to determine the long-term behavior of the exact 
solution. 

Medium Forcing 

For all the medium-forcing cases, as for all the large-forcing cases, the nonphysical long-term motion 
was assumed at large times. However, the short-term and transient solutions of the medium-forcing cases 
differed from those of the large-forcing cases, and they persisted for much greater times. For example, 
consider the medium-forcing, low natural-frequency case (« = 0.2033 k a ) plotted in figure 10. Initially, 
the motion is a simple, harmonic motion with a frequency equal to the Strouhal frequency. As the motion 
slowly decays, the characteristics of the long-term solution slowly appear: the frequency of the motion 
remains equal to the Strouhal frequency, odd harmonics of the Strouhal frequency contribute more to 
the solution (fig. 10b), the velocity gradually develops jumps, and the acceleration and forcing function 
gradually develop pronounced peaks at intervals 8t = n/k a . The peak amplitude of the acceleration is 
equal to the amplitude of the forcing function; this result agrees with that of the long term solution. At no 
time are high-frequency oscillations observed in the computed solution. 

The position and velocity, however, do not show the f -1 / 2 decay rate characteristic of the long- 
term motion (fig. 10c); instead the decay rate appears to vary with time. For instance, at the beginning of 
the computation, the motion does not appear to show any signs of decay; at the end of the computation, 
the amplitudes of the position and velocity are proportional to t~° 4 . As the computation proceeds and 
the jumps in the velocity become more pronounced, the solution will better approximate the conjectured 
long-term solution and the position and velocity will be proportional to f -1 / 2 . 

The smooth, continuous transition to the long-term behavior is reflected in the smooth behavior of 
the Lyapunov exponents (fig. lOd). All three exponents are negative implying that the phase space defined 
by ( h, h, k a tV\ + h 2 ) is contracting. The two small Lyapunov exponents appear to be approaching zero. 
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and, again, this may be an artifact of the slow convergence rate of the solution as opposed to the solution 
lying on a two torus. At large times, the large, negative Lyapunov exponent and the sum of the exponents 
is continuously decreasing at a rate proportional to Kf^ = —0 .001 1° 69 . The sum has not yet assumed the 
t x / 2 behavior of the long-term solution. This is not unexpected since the velocity has also not assumed its 
long-term behavior. 

At any time the sum of the Lyapunov exponents for the medium-forcing case is much smaller 
than the sums for the large-forcing cases. The dominating term of the sum at large times is proportional 

to theCL cos (^k 3 tV\ + A 2 j . Both h and eCx are smaller for the medium forcing case; therefore, the sum 
of the Lyapunov exponents should be smaller. As t — * oo and the computed motion better approximates 
the conjectured motion, h for the medium-forcing case will be identical to h for the large-forcing cases, 
the sum of the Lyapunov exponents will be proportional to t 1 / 2 , and the ratio of sums of the Lyapunov 
exponents for the medium- and large-forcing cases will be equal to the ratio of the forcing amplitudes. 

For the medium-forcing and matched-frequency case, the transition to the long-term solution is 
smooth and has many of the features just discussed for the medium-forcing, low natural-frequency case; 
the velocity gradually develops jumps and the amplitude of the motion is proportional to t~P (figs. 11a and 
lib). Since h for the matched frequency solution is generally larger than for the low frequency solution, 
there is a small modification to the acceleration and forcing function. The larger A causes additional wig- 
gles in the forcing function and acceleration, as is observed in figure 11a. These plots, however, do show 
the same development of pronounced peaks and rapid changes at 8t = n/k a intervals and the additional 
wiggles in the solution slowly die out. For this case, the second Lyapunov exponent does not remain small 
(fig. 11c). However, at the end of the computation, it does become smaller and appears to be asymptoting 
to 0. 


When the natural frequency is slightly larger than the Strouhal frequency, the transition to the long- 
term behavior shows features other than those discussed previously. The Fourier transform of the initial 
motion has peaks at the Strouhal frequency and its odd multiples, at the natural frequency and its odd 
multiples, and at side bands with a spacing equal to the difference in the Strouhal and natural frequency. 
These multiple frequencies result in beating in the position and velocity plots (figs. 12a and 12b). 

As the time increases, the amplitude of the motion decreases, Fourier terms with the Strouhal fre- 
quency and its odd multiples begin to dominate the motion, and the beating disappears (figs. 12b, 12c, 
and 12d). However, the motion still has not yet made the transition to the long-term behavior. Instead of 
minimization of the acceleration and forcing terms, there are high-frequency oscillations in both terms. 
The frequency of these oscillations increase until the system makes some attempt to minimize the accel- 
eration and forcing function. At this point, the motion does decay to the long-term behavior in a manner 
similar to that of the large-forcing, large natural-frequency case. 

The first and third Lyapunov exponents for the medium-forcing, large natural-frequency case are 
well-behaved; the first is small, positive, and approaches 0 , the thirdis large, negative, and has a t? behavior 
(fig. 12e). The second exponent is negative and remains negative; however, there is a sudden change in its 
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slope. This may be related to the loss of beating and the transition to a motion that is composed of Fourier 
terms with frequencies equal to the Strouhal frequency and its odd multiples. 

The short-term motions for these three medium-forcing cases do show some of the features of ex- 
perimental data; the motions are periodic and beating is observed, at least initially, for one of the cases. 
However, they do not adequately describe the experimental motion. For example, in all three cases, a short- 
term motion exists that is periodic, has a frequency equal to the Strouhal frequency, and does not minimize 
the forcing function; however, Feng’s (ref. 4) data seems to indicate that the frequency of the motion and 
the shedding frequency should not be locked together for at least one of these cases, the k = 0 .2033 k s 
case. Also, for all three cases, the short-term calculated motion slowly decays; the experimental motion 
does not. Hence, comparisons between calculated and experimental amplitudes are impossible. Further- 
more, the initial motion for the k = 1 . 1 k a case shows signs of beating; Feng (ref. 4) makes no mention of 
beating for these parameters. 

All three medium-forcing cases eventually do assume the nonphysical long-term behavior. The 
length of time required before the motion settles into the long-term behavior is much larger than that for the 
large-forcing cases. This results directly from the equation of motion’s inability to support large magnitude 
motions; for the large-forcing case, the initially large h causes high-frequency oscillations in the forcing 
function and the motion damps out. For the medium-forcing case, the amplitude of the motion is much 
smaller, the forcing function is not forced through many oscillations, and the transient can maintain itself 
over an extended period of time. 


Low Forcing 

For the low-forcing amplitude, two computations were performed, one with the natural frequency of 
the spring equal to the Strouhal frequency and the other with the natural frequency equal to 1.1 times the 
Strouhal frequency. Neither of these cases reached the long-term behavior; this is undoubtedly due to the 
small amplitudes of the motions. As observed in the large- and medium-forcing cases, a decrease in the 
forcing amplitude can greatly delay transition to the long-term behavior. 

For the matched frequency case, the motion is damped at a t -1 / 2 rate and the Fourier transforms 
show peaks at the Strouhal frequency and its odd multiples (fig. 13). This behavior is in agreement with 
the damping rate and Fourier transforms of the long-term motion. However, unlike the long-term motion, 
the subharmonics of the Strouhal frequency for the low-forcing, matched-frequency case have a negligible 
effect on the motion, and the plots of the position, velocity, and acceleration show that the motion is 
sinusoidal and that there is no minimization of the forcing term or acceleration. Since the ratio of the 
amplitudes of the higher order harmonics to that of the fundamental term do not change appreciably with 
time, a possible solution appears to be one of a slowly damped, approximately sinusoidal motion. 

The Lyapunov exponents for this low- forcing, matched-frequency case are plotted in figure 13d. 
There are two small exponents that do not appear to be asymptoting to any particular value. Not only is 
there a lot of variation in these exponents, there are sudden changes in their general direction. The sum of 
these two small exponents does appear to be asymptoting to a negative value. 
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The largest negative exponent and the sum of the exponents are not as large as those of the large- 
forcing and medium-forcing cases. They do, however, share some important features, namely, the sum of 
the exponents is dependent upon time and the size of the sum is dependent upon eCi- For this particular 
case, the computed sum of the Lyapunov exponents is given by 2 .8 x 10 _5 i° 54 . The exponential behavior 
of the sum is in good agreement with the theoretical value of the sum for the long-term behavior, and the 
coefficient is of the order of, but not equal to, the coefficient predicted for the long-term motion. There 
should be some difference in these two values since the computed motion is sinusoidal and has not assumed 
the long-term behavior. 


For a sinusoidal motion, A = At~ x ^ sin( k a t + <p ) , the sum of the Lyapunov exponents can be written 

as 

/ k A 2 \ 

^ A, = eChk s At' /2 sin(k a t + <p) cos i k s t + — cos(2k a + 2 <f>) J 

where only the dominant term of the trace of the Jacobian is considered and A 2 has been assumed to be 
much less than one. Taking the average on the interval to — it jk a < t < to + n/k a gives 


X>= j*C L k a At l/2 


A 2 


cos ( (f>) Jo ( k 3 — ) + sin(</>) J] ( fc a — J j 


A 2 Y 


where the J t ’s are Bessel functions. For this computation, tCi = 0 .00007, A & 5 .37 , cos (</>) ^ 1/3, 
and the analytic sum is 


= (0.2 x 10~ 5 cos(</>) +4.4 x 10" 5 sin (</>)) t 1/2 
ps 4 .0 x 10- 5 t 1/2 

If the analytic and computed sums are to be compared, this expression must be divided by 1 .5 (see section 
discussing computation of the Lyapunov exponents), or 

y^A, = 2.7 x 10- 5 t 1/2 

This compares favorably with the computed sum; differences between the two sums may be due to the 
approximations of the amplitude of the motion and its phase angle and by the approximation of the trace. 


In an attempt to understand this damped, very long-lived transient and to determine whether this 
motion can persist indefinitely, assume that the motion has a small amplitude and is periodic, h = 6 sin ( tot + 
(p) . Substituting into the forcing function and assuming that < 1 gives 


sin (^k a t\/ 1 + A 2 j = sin A: s f\/l + w 2 6 2 cos 2 (u>t + <£)j 

/ g\2 

k a t + k a t f y- j cos(2wf + 2<p) 


= sin 
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The Fourier transforms of the forcing function are written as 


F .T .(cos) = 


t\ — to 


k a t + k a t l— J cos(2wf + 2<£) cos(kt)dt 


F .T .(sin) = 


2 r ti . 

= / sin i 

1 1 — *o Jto 


k a t+k a tl—J cos(2wt + 2(f>) sin (kt)dt 


2 

If the integration interval t\ - to < fo, then k a t (Jj-) approximately constant over the integration 
interval, and the Fourier transforms can be written as 


F.T.(cos) = 


ti — to 


sin [k a t + /?cos(2k>£ + 2 (/>)] cos(kt)dt 


F.T.(sin) = 


ft — fo 


sin [k a t + /?cos(2u;f + 2 <f>)] sin(kt)dt 


where /? = k a to(u8/2) 2 . Although h is small, (d is not small if fo is large, and the /? cos ( 2 u;t + 2<f>) term 
cannot be ignored. Letting r = 2wf + 2<f> gives 


F.T.(cos) = 


ri — 70 


-(t - 26) + /3cos(t) cos - — (T—2<f>) dr 
' 2u 


T\ ~ To 


(r- 2<f>) + (3cos(t) sin - — (r — 2<f>) \ dr 


By assuming k a /u> = N/M where N and M are integers, it is possible to determine the Fourier transforms 
analytically. The forcing function will repeat itself after an interval A r = 4 M n, and, therefore, the 
Fourier transform will be zero unless k/2u> takes on values 1/2 M, 2/2 M, . . . , P/2 M, etc. Letting 
ti = 70 + 4 M 7t and integrating gives 

For 2ft 

F.T.(cos) = 0 
F.T.(sin) = 0 

For ^ = 2ft, 2ft 

FT. (cos) = sin -2 ft/) J a (/?) 

-2 ft*) Ja, ()9) 
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¥or *L=JL ilRu U+IL = 2 R 2 


F.T.(cos) = sin — -2R 2 <^j J Rz (P) 

F.T.(sin) = - cos ^^-2 R 2 ^ J Rl (/ 3 ) 

Yox^ = 2Ru^f = 2Ri 

F.T.(cos) = sin — ^R\<^j Jit, (P) + sin - 2 # 2 *^ Jit 2 (P) 

FT. (sin) = cos (^-2R ]( ^j J fll (/3)-cos (^-2 R 2 *) J r 2 (P) 

where P is an arbitrary integer and the J^’s are Bessel functions of order R. The amplitudes of the various 
Fourier terms are 

For B&12Ru*ff-12R 1 

Ap = 0 

Fot^ = 2R u ^¥2R 2 

Ap = Ja C/5) 

For^£^2ft. i % £ = 2H 2 

Ap = J r 2 (P) 

For ^fr~ = 2 Ru - 2 R 2 

1 /2 

Ap = I Jjjj (P) + J^ 2 (P) —2 cos 

Consider the case where u is equal to the Strouhal frequency, k s . The only terms of the Fourier 
expansion of the forcing function that will have nonzero coefficients are those that satisfy P/M = k/u = 
2 R + 1 where R is any integer. Therefore, the Fourier transform of the forcing function only has terms 
with frequencies that are odd multiples of the Strouhal frequency. The amplitudes of these terms are 

A r = [j l R {p) + J| +1 ( P ) - 2 cos (2<t> - 0 J - R (P) Jjr + i ( P )\' 2 


20 



If / 3 is large, this expression can be written as 


Ar =J— [l + sin(20) cos(2/?)] 1/2 


nonlinear equation into a linear equation by assuming that h <C 1 gives, for u = 

OO 

(2 ij, — cCd ) h+ kfh = cCl ^ Ar sin [ ( 2 R + 1) k a t + 6 r] 

r = o 

where the Ar ’ s are the amplitudes obtained from the Fourier transform of the forcing function and the 6r s 
are the phases for these terms. Expanding the motion into a Fourier series, 

OO 

h = ^Br sin [(2 R+ 1 )k a t + (j> R ] 

R = o 


Transforming the 

K kgy 

h + 


and using the linearized equation of motion to solve for the coefficients Br gives 

[l + sin(2</>)cos(2/3)] 1/2 4eC L 


Br = 


k aS /(JW- + 4 RyHZ + (2/i - cC D ) 2 (2R + l) 2 


1 P =: 

V 2i[k a to 


The ratio of Br/ Bo is 


Br (2/x — cCp) 

Bo " ^(4^ + (2/i-cCD) 2 (2fl+ l) 2 

/2/i — cCd\ 1 

) R 2 + R 

The ratio of the amplitudes of the harmonic terms is independent of time, and, since the total damping is 
small (2 .2 x 10 -4 ), the higher-order harmonics are negligible when compared to the fundamental term, 



4.4 x 10” 5 


1 

WTr 


The assumption that the motion is sinusoidal is therefore reasonable. 


The analytic results just discussed are in good agreement with the computed results. Both the anal- 
ysis and computation show that the motion is primarily sinusoidal with a frequency equal to the Strouhal 
frequency. The Fourier transforms of the analytical and the computational results do show peaks at odd 
multiples of the Strouhal frequency, but these have amplitudes several orders of magnitude less than that 
of the Strouhal frequency. The analytic and computational ratios of the amplitudes of the third and fifth 
harmonic with respect to the amplitude of the fundamental are in good agreement. There is some difference 
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in the analytic and computational ratios of the higher order harmonics; this may be due to numerical limita- 
tions, errors caused by linearization of the equations of motion, deviations from the large /5 approximation, 
or small contributions to the forcing function from the higher-frequency components of the motion. 


The analytic amplitude of the fundamental term, 8 = Bo, can be written as 


Bo = 


1 eC L 

"1 + 2 sin(2</>) cos(2k a (k 3 Bo/2) 2 to)~ 

/ 2/4 - cCd 

2 'nk a to 


1/4 


In figure 13c, the amplitude of the motion is plotted versus time. Note that the motion does die off at a rate 
r^ - 1 / 2 . This implies that /5 = k s to ( k s Bo /2) 2 is constant and that the phase of the motion has to vary in 
a manner such that 

{1 + 2 sin (2(f>) cos[2k 3 (k a Bo/2) 2 t 0 ]y /4 = 


or 


sin(2<£) = 


1 


cos (2/5) 


t 


- 1 


The computations did show very minor variations in the phase angle; hence, the analytic and computed 
results are not in disagreement. 


These analytic and computed results do share one similarity with Van Atta and Gharib’s (ref. 2) 
experiments. In their experiments, a wire with one of its subharmonics tuned to the Strouhal frequency os- 
cillated with a frequency equal to the Strouhal frequency. Unfortunately, their Fourier transforms were not 
taken over a wide enough range of frequencies to capture any higher order harmonics of the Strouhal fre- 
quency. The experiment and computation do show disagreement concerning the amplitude of the motion; 
the computed motion showed signs of decay, the experimental motion did not. 

There are no indications in either the analysis or in the computations that the damped, nearly sinu- 
soidal motion will give way to the nonphysical, long-term motion. This slowly damped, nearly sinusoidal 
motion, however, can not persist indefinitely. As time increases, the amplitudes of h, h, and h will continue 
to decrease until they are of the order of the forcing term. For the low-forcing, matched-frequency case 
discussed here, this will not occur until t m 10 8 — 10 9 . At this point, a nearly sinusoidal motion will 
not satisfy the equations of motion and the solution will have to assume another behavior, presumably the 
long-term behavior, as both the large- and medium-forcing cases did. 

For the mismatched frequencies, the above analysis does not hold. If the motion oscillates primarily 
at the natural frequency k (as indicated by the computations, fig. 14a), the forcing function will have Fourier 
terms with frequencies k = 2Rn - k s and 2 Rk + k s , R is an arbitrary integer. None of these terms will 
force the motion at a frequency k unless k can be written as k = k a /R; hence, a simple harmonic motion 
cannot be assumed. If, on the other hand, we assume that the motion is a two frequency motion and can 
be written as 

h = At -1 sin ( k a t + fa ) + Bt~P cos ( nt + fa) 
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then for large times and short time intervals, the forcing function can be approximated by 


sin ^A; a t\/l + 2 j = sin{fc s f + acos(2k a t + 2 <f>\) + /?cos(2/cf + 2 <f> 2 ) 

+ 7COS [(/c s + n)t + (f) 1 + </> 2 ] + 6 cos [(/c 3 — «)£ + </»i— </>2 ] } 


(- 1 )’ 


= Y' J 4„sin((2( + 1 )h,t + 6i) Y) ' —- {ffcosflKt + 2^2 ) 

M £S( 2m)! 


+ 7 cos [(k a + n)t + <f>\ + <f >2 ] + 8 cos [ ( — k) t + 4>\ — <f >2 ] } 


2m 


E®. sin[( 2 /+ 1 )k 9 t + e' l \ 


(- 1 )’ 


i=0 


m=0 


(2m + 1) ! 


[/5cos(2 Kt + 2<f>2) 


+ 7 cos [(k s + n)t + (f>\ + <f> 2 ] + 8 cos [( k s — /c)f + </>i — 0 2 ] } 


2 m+l 


The sum over m introduces frequencies that can be written as 2 Ln + M( A: s — «) + iV( /c 3 + «) where L, 
M, and TV" are integers, and the frequencies produced by the forcing function are given by (2 K + M + 
N+ l)k a + (2L + N-M)K = (2K+2L+2N+1 )k+( 2K+ M + N+ l)(k a -K). This implies that 
the Fourier transform of the motion will have terms at the natural frequency, its odd multiples, and at side 
bands with the spacing given by k — k a . It is impossible to obtain the relative amplitudes of these terms 
through simple analysis; however, the computed Fourier transforms for the « = 1.1 k a case do show peaks 
at the natural frequency, its odd multiples, and at side bands (fig. 14b). Not all of the analytically allowed 
frequencies are observed in the computations; their amplitudes may be so small that they are lost in the 
background noise. 


The peak amplitudes of the computed position and velocity for the low-forcing, unmatched- 
frequency case (fig. 14c) did not decrease monotonically as did the amplitudes for the matched-frequency 
case. This is undoubtedly due to the interaction between the multiple frequency terms and time in the 
argument of the forcing term. Although complete analysis of the Fourier expansion of the forcing func- 
tion is too difficult to determine the peak amplitude of the motion as a function of time, it is expected that 
the amplitude, on the average, will continue to decrease, and at some point this case will too go into the 
long-term behavior. 


These results do share some similarities with Van Atta and Gharib’s data (ref. 2), although there 
are some discrepancies. For mismatched natural and Strouhal frequencies, the experimental data did show 
peaks at the natural frequency, the Strouhal frequency, and at sidebands. However, in the experimental data, 
the peak at the Strouhal frequency dominated; in the computed results, the peak at the natural frequency 
of the spring dominated. Furthermore, the experimental data showed no decay in the amplitude, whereas 
the computed amplitude did slowly decay. 


The Lyapunov exponents for this case are all small (fig. 14d). The first exponent is small and positive 
and shows no indications of becoming negative; it does, however, appear to have a zero asymptote. The 
presence of this small positive exponent is indicative of chaos. Although the motion does not appear 
particularly chaotic, changes in the amplitude and spectra of the motion may be chaotic; this could result 
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in a positive Lyapunov exponent. The second exponent is small, negative, and also appears to have a zero 
asymptote. The third exponent is also small, negative, and will remain negative. 

Contrary to the results of the previous cases, there is no large, negative Lyapunov exponent, and the 
sum of the exponents do not have the characteristic t 1 / 2 behavior. Unlike the previously discussed cases, 
the frequency of the motion has not locked onto the Strouhal frequency and the amplitude of the motion 
is not dying off at the f -1 / 2 rate. Therefore, the approximations used to analytically verify the behavior 
of the sum of the Lyapunov exponents are not valid here, and the characteristic t 1 / 2 behavior will not be 
observed until the motion assumes the long-term behavior. 

CONCLUSIONS 


The nonlinear, quasi-static model developed by Tobak et al. (ref. 5) to describe the motion induced 
by vortex shedding of a cylinder has some interesting features, although it does not adequately model the 
problem. The computed short-term motion (which may persist for thousands or millions of periods) shows 
several of the features of the experimental motions including the locking of the frequency of the motion 
to the shedding frequency, periodic motions, and beating. Unfortunately, the parameters at which lock-in, 
periodic motions, and beating occur do not always agree with experiment, and the weak damping of the 
short-term computed motions is not observed experimentally. 

Computations and analysis indicate that the long-term motion is independent of the lift, the drag, 
and the natural frequency of the cylinder. This long-term motion is periodic with a frequency equal to the 
Strouhal frequency (the shedding frequency), and it minimizes the forcing function. As a result of this 
minimization, the acceleration is nonzero only at pronounced peaks of finite height that decrease in width 
as time increases, the velocity exhibits rapid jumps, and the motion of the cylinder decays with a t -1 / 2 
rate. This weakly damped, long-term motion does not agree with experimental results that indicate that 
the motion persists and that the amplitude of the motion is dependent upon the Reynolds number and the 
natural frequency of the cylinder. 

This nonphysical motion is directly attributable to a motion dependent forcing frequency that intro- 
duces a dependency on the time origin into the equations of motion. It also effects the Lyapunov exponents 
and results in their sum and at least one of the exponents being unbounded. Furthermore, the slow damp- 
ing rate of the motion and its strange long-term behavior make it impossible to determine the character 
of the asymptotic motion from the Lyapunov exponents. For all cases, the long-term motion involves a 
periodic switching back and forth between two fixed points. The computed Lyapunov exponents for the 
cases studied indicate that the long-term behavior may be a fixed point, lie on a one- torus, or even lie on a 
two-torus; hence, they cannot be used to adequately define the long-term motion. 

This one-degree-of-freedom model does not describe the motion induced by the vortex shedding of a 
cylinder; adding a time-dependent phase term or a second degree of freedom may prevent the minimization 
of the forcing function. Then, a motion that did not settle into the nonphysical behavior of this one-degree- 
of-freedom system might be obtained, as well as more favorable comparison to experiment. 
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Figure 3 - Large forcing amplitude, small natural frequency, (a) Motion. 
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Figure 3 - Concluded, (b) Phase portraits, (c) Amplitude of the motion: h- \ .19t 0 5 , /i = 2 .19 1 
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Figure 4 - Conjectured velocity. 
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Figure 5 - Lyapunov exponents for the large forcing amplitude, small natural frequency. 
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Figure 6.- Concluded, (b) Amplitude of the motion, (c) Lyapunov exponents 
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Figure 7 - Variation of step size and number of significant digits. 
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Figure 9 - Effect of starting amplitude on position, velocity, and acceleration. 
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Figure 10- Medium forcing amplitude, small natural frequency, (a) Motion. 
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Figure 10.- Concluded, (b) Fourier transforms, (c) Amplitude of the motion. For the last three points 
h = 037t~ OA /h = 0 A9t~° 4 . (d) Lyapunov exponents. 










Figure 11- Medium forcing amplitude, matched frequencies, (a) Motion. 
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Figure 11.- Concluded, (b) Amplitude of the motion, (c) Lyapunov exponents. 
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Figure 12.— Medium forcing, k = 1 .1 k 3 . (a) Short-term motion, (b) Fourier transforms. 
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Figure 12.- Concluded, (d) Amplitude of the motion, (e) Lyapunov exponents. 
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Figure 13 - Concluded, (b) Fourier transforms, (c) Amplitude of the motion: h = 4 
h = 5 ,37 -0 51 . (d) Lyapunov exponents. 
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Figure 14 - Small forcing amplitude, rz = 1 . 1 k s . (a) Motion. 
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